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We introduce a versatile method to compute electronic steady state properties of strongly corre- 
lated extended quantum systems out of equilibrium. The approach is based on dynamical mean-field 
theory (DMFT), in which the original system is mapped onto an auxiliary non-equilibrium impurity 
problem imbedded in a Markovian environment. The steady state Green's function of the auxil- 
iary system is solved by full diagonalization of the corresponding Lindblad equation. The approach 
can be regarded as the nontrivial extension of the exact-diagonalization based DMFT to the non- 
equilibrium case. As a first application, we consider an interacting Hubbard layer attached to two 
metallic leads and present results for the steady-state current and the non-equilibrium density of 
states. 



Due to the progress made in microscopically control- 
ling quantum mechanical systems within quantum op- 
tics and ultracold quantum gases drS], in solid state 
nanoscience, spintronics, molecular electronics, [6] as 
well as ultrafast laser spectroscopy [71-fTO]. the interest 
in correlated systems out of equilibrium has steadily 
increased in recent years. These achievements have 
prompted new and boosted old related theoretical ques- 
tions such as nonequilibrium quantum phase transi- 
tions pj] , dissipation and decoherence [12] , and thermal- 
ization after a quantum quench [T5l415| . 

In this respect, the theoretical description and under- 
standing of strongly correlated quantum systems out of 
equilibrium poses an exciting challenge to modern the- 
oretical physics. A widely used and successful method 
to treat strongly correlated lattice systems in equilib- 
rium is dynamical mean-field-theory [I5HIS] (DMFT). 
The success of the method lies on the one hand in the 
nontrivial treatment of dynamical properties, and on the 
other hand in its applicability to a range of different 
problems, from solid states fermionic systems to ultra- 
cold bosonic atoms, as well as the possibility to com- 
bine it with realistic electronic structure methods [!§] . 
Recently, DMFT has been extended to deal with time 
dependent nonequilibrium problems [20-25]. The exten- 
sions are based on the Kadanoff-Baym-Kcldysh nonequi- 
librium Green's function approach [26lf29] . 

DMFT relies on the solution of a correlated impu- 
rity problem, which constitutes the bottleneck of the 
approach. Several techniques have been adopted in the 
equilibrium case. Most of them have been applied, in a 
more or less approximate or limited way, to nonequilib- 
rium DMFT as well, either in steady state or within full 
time dependence. Among them are iterated perturbation 
theory |20) . numerical renormalization group [23 , con- 
tinuous time quantum monte carlo (CTQMC) [2~4"l 150], 
noncrossing approximation and beyond |31j . Addition- 
ally, exact DMFT solutions are available in certain lim- 
its |21|, I32[ I33j . Nonequilibrium quantum impurity 
problems (not within DMFT) have also been studied 



by means of scattering-states approaches [3U |35J, per- 
turbative methods [361 137] in combination with renor- 
malization group (RG) [3H , time-dependent density- 
matrix RG [40l Hi] and numerical RG [42] flow equa- 
tion [43], functional RG (44] [45], dual fermions [46], and 
finally CTQMC on an auxiliary system with an imagi- 
nary bias [47PS] , 

In this Letter, we propose an approach that, in con- 
trast to previous work, while directly accessing steady- 
state properties, features a solution of the DMFT im- 
purity problem with controlled accuracy. This means 
that the accuracy can be directly estimated by comparing 
the exact and the effective bath hybridization functions 
(Fig. [2j. Also, no often unreliable analytical continua- 
tion is required. At the heart of the method lies a 
solution of the nonequilibrium DMFT impurity problem, 
which can be seen as a generalization of the exact diago- 
nalization (ED) approach, widely used in the equilibrium 
case |16) . However, a crucial difference with respect to 
conventional DMFT-ED is the fact that here the effective 
impurity model describes an infinite system and, thus, 
displays a continuous spectrum. 

In ED-based equilibrium DMFT [TB] a certain number 
of noninteracting bath sites is introduced in order to fit 
the bath hybridization function required by the the self- 
consistency condition. The maximum number of bath 
sites is limited by the exponential increase of the many- 
body Hilbert space. In equilibrium, the fit is carried out 
in imaginary (Matsubara) frequency space, where func- 
tions are smooth, in contrast to real frequency. There are 
several difficulties in this approach when trying to extend 
it to nonequilibrium steady states, (i) Due to the finite 
number of bath sites, a stationary solution of the impu- 
rity problem will always produce some equilibrium self 
energy. Besides the fact that this may be questionable, 
it is not clear which chemical potential or temperature 
should be used for the impurity problem, (ii) Secondly, 
due to the finite system, the bath spectrum is discontinu- 
ous, so that a fit in real frequencies becomes problematic. 
Unfortunately, there is no Matsubara Green's function in 
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Figure 1: (Color online) Schematic representation of the sys- 
tem at study. 

nonequilibrium, so that this poses a serious problem. 

The alternative presented in this Letter consists in re- 
placing the DMFT impurity Hamiltonian with an effec- 
tive one which is solvable by ED but at the same time 
describes a truly infinite system. This is obtained by con- 
necting the interacting impurity to a moderate number 
of bath sites which, in turn, are attached to Markovian 
reservoirs, see below for details. The exact bath spec- 
tral function is smoothly obtained in the (ideal) limit of 
an infinite number of bath sites. The action of such 
Markovian baths on the reduced density matrix of the 
system consisting of the other bath sites and of the impu- 
rity, is described by the Lindblad quantum Master equa- 
tion [50] . The latter can be readily solved by diagonaliz- 
ing the Lindbladian within the many-body "super-Fock" 
space of reduced density matrices of the system. Its solu- 
tion determines both the retarded and Keldysh impurity 
Green's function, as well as the self-energy. The latter is 
used in the DMFT loop to obtain the new bath hybridiza- 
tion function, which is fitted by new bath parameters. 

In order to illustrate the approach, we apply it to a 
simple model describing a heterojunction consisting of a 
correlated interface (c) sandwiched between two metal- 
lic leads a = l,r (see Fig. [I]). Experimentally, such a 
setup has been recently explored where the correlated 
layer was realized by a V2O3 microfilm that is coupled 
to Au leads [5T] . 

Before a certain time r < r the three regions c, I, r are 
disconnected and in equilibrium at different chemical po- 
tentials Hct^Ii and /i r , respectively. This amounts to ap- 
plying a bias voltage $ = \i\ — fi r between the leads. The 
central region, lying on the x — y plane, is described by 
a single-band Hubbard layer on a square lattice with on- 
site interaction U, onsite energy e c , and nearest-neighbor 
hopping t. The leads consist of two half-infinite cubic 
lattices described by a nearest-neighbor noninteracting 
tight-binding model with hopping which we take as unit 
of the energy t = 1, and on-site energies e a . We re- 
strict for simplicity to the particle-hole symmetric case 
for which e c = —U/2, fi r = —fJ,i, and e r — —£i- Fi- 
nally, we take /i a — s a , which corresponds to having 



the same electron densities in the two leads. A related 
nonequilibrium model has been treated in DMFT within 
perturbative impurity solvers in Refs. |25[ 131] . 

Starting at t = t , a nearest-neighbor hopping v is 
switched on between the central region and the leads. Af- 
ter a sufficiently long time, a steady-state is reached, pro- 
vided no trapped surface states occur. Nonequilibrium 
properties, in general, and nonlinear transport in partic- 
ular can quite generally be addressed in the framework 
of the Keldysh Green's function approach [26 29, 52 1. 
Here, we adopt the notation (see e.g. Ref. j^S]) where 
the (underlined) Keldysh Green's function is a 2 x 2 ma- 
trix containing the retarded (G R ), advanced (G A ), and 
Keldysh (G ) components. In steady state, these depend 
on a single frequency only. The system is translation in- 
variant in the direction parallel to the layer, which we 
denote as ||. Accordingly we can write Dyson's equation 
for the layer (c) Green's function G(ky,£j), ky being the 
|| momentum [52 as 

am,^- 1 = g-\u)- ]r « 2 £Q (k||,c)-s(k||, w ) . (i) 

Here, £(ky,u;) is the self-energy, g Q is the v = 0, U = 
layer Green's function, and g a (k.^,u>) are the v — 
leads Green's functions on the first lead layers. Their 
retarded and Keldysh components are readily obtained 
analytically in terms of the Green's function of a half- 
infinite tight-binding chain. 

Within DMFT, one approximates the self energy by 
a local, i.e. ky -independent S(w), which is determined 
by solving a quantum impurity model with the same in- 
teraction U embedded in a self-consistently determined 
bath [IB]. The latter is completely specified by the bath 
hybridization function A(w), which is determined self- 
consistently by requiring that the Green's function of 
the impurity G IMP (u;) = (Qoiuj)^ 1 — £(u;)) be equal 
to the local Green's function of the layer (cf. [TBI |2"T1 |2"S"] ) 
GlocM = I (Sja G(k||,o;), where G(ky,w) is given 
by ([!]) with E(ky,tL>) = S(w) as obtained by the solu- 
tion of the impurity problem. Here, Go is the "Weiss" 
bare Green's function of the impurity model defined as 

SoM" 1 ^^)- 1 -£(«). 

The solution of the impurity problem is in fact the 
DMFT bottleneck. The usual (equilibrium) DMFT ED 
procedure consists in approximating the effect of the to- 
tal bath hybridization function A(w) by an "effective" 
bath with a finite number TVf, of bath sites. Quite gen- 
erally, in equilibrium one carries out some fit to the 
bath hybridization function in Matsubara space. As 
discussed above, this is not appropriate in nonequilib- 
rium steady state. In this Letter, we present a differ- 
ent approach: In additional to a certain (even) number 
Nb of bath sites, which we more conveniently connect 
to the impurity in the form of two chain segments, we 
include two Markovian baths, which represent a par- 
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ticle reservoir and sink, respectively. Their role is to 
compensate for the "missing" part of the infinite chain 
which would be necessary to exactly reproduce the de- 
sired A(w). The bath parameters, i.e. hopping and 
on-site energies of the bath sites, as well as the Lind- 
blad coefficients (see below) of the Markovian baths, are 
then fitted to A(cj). More specifically, we minimize the 



cost function J dio (Im(A a: (a;) 

x=R,K 

A(w) is obtained from G hoc via 



A effM)) , where 



AM = ^H- 1 - GlocH- 1 - E(w) 



(2) 



while A cff is the bath hybridization function produced 
by the effective bath (bath sites+Markovian baths). An 
important aspect is that, although the (outermost) baths 
are Markovian, their effect on the impurity site is non 
Markovian due to the presence of the intermediate bath 
sites. This can be seen, for example, in the spectrum of 
A off in Fig. [2j which in the Markovian case would be a 
constant. Furthermore, upon increasing the number Nj, 
of intermediate bath sites, the effect of the Markovian 
bath becomes weaker and one is expected to approach 
the exact result A eff (w) — A(w) for large Nb- 

We now specify the effective bath more in detail. This 
consists of an Hamiltonian for the "system" (a chain of 
impurity+bath sites) 



H 



E R - 



^ C 0T C 0T C (H C 0i> 



(3) 



in usual notation, where is the impurity site, and 
n = —1, • • • , — I (n = 1, • • • , I) are left (right) bath sites. 
Here, Eq^ = e c , and the bath energies E n>n , n ^ as 
well as the hoppings E„ tm , n ^ m are fit parameters |53j . 
whereby one can restrict to nearest-neighbors E nn ±i. 
The effect of the Markovian baths is expressed in terms 
of the Lindblad quantum master equation which controls 
the time (r) dependence of the reduced density matrix p 
of the system [50]: j^p — Cp, where C = Ch + £&> an d 
ChP = — i[H,p\. The dissipator Cb has the form 

C-bP = 2 y ] (j^n]m { c no-P c tncr ~ C " lo " C ™°"^ ) 

+r (2) pc --\oc J \Y\ 

' x n.m K^nar^ma nlri'-'lTUT naj J J ' 

with real, symmetric, and positive definite Lindblad ma- 
trices r^ m . At first sight, one would assume a "source" 
bath attached to the leftmost and a "sink" bath to the 

rightmost site. However, the accuracy improves consid- 

(i) 

erably if one allows all T n 'm to be nonzero parameters to 
fit A. 

As discussed in detail in Refs. [54 57 , the open- 
system problem describing the effective bath can be 
mapped onto a superhamiltonian acting on a super- 
fermion space with twice as many degrees of freedom (i.e., 
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Figure 2: (Color online) Imaginary part of the DMFT 
and effective bath hybridization functions (retarded (R) and 
Keldysh (K) components) for $ = 2, v 2 =0.1, and two differ- 
ent values of U. We also plot the effective bath hybridization 
function (A cff 4 ) obtained with Nt = 4. 




Figure 3: (Color online) Imaginary part of the impurity and 
local Green's function for (7 = 1 (left) and (7 = 4 (right) and 
$ = 2. 



"orbitals"). This many-body superhamiltonian, corre- 
sponding to iC, which is non-hermitian, can be diagonal- 
ized by conventional methods within the super-Hilbert 
space. Quite generally, L has a unique eigenvector with 
eigenvalue 0, which corresponds to the steady-state den- 
sity matrix pss- All other eigenvalues have a nega- 
tive real part, corresponding to decaying terms. With 
the same formalism, and exploiting the quantum regres- 
sion theorem |58) . one can evaluate correlation functions 
Cab(t) = trA(T)B(0)pss of any pair of system opera- 
tors A and B, and thus the required impurity self energy 
E(w). [59] 

The noninteracting Green's function for the effective 
system+bath is necessary in order to extract S(w) and 
the bath hybridization function A cff . This can be easily 
obtained [50] by observing that the Markovian baths can 
be exactly represented by two noninteracting fermionic 
baths in the wide-band limit with chemical potentials 
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±oo. By taking into account the relation between the ma- 
trices r and the parameters of this bath [61] , one obtains 
for the noninteracting system Green's function (boldface 
object represent matrices with indices corresponding to 
system sites n): {G^) R = ul - E + i(r« + r( 2 )), 

(G- 1 ) K = 2i(r( 1 )-r( 2 )). 



The DMFT self-consistency loop consists in (i) start- 
ing from some initial values of the variational parameters 
E n>m and Tn}m, (ii) solving the impurity problem via the 
approach described above and determining E, (iii) eval- 
uating G LOC (u}) and A(w) (from Rib) (iv) determining 

(i) 

new values of the parameters E n ^ m and r„; m by mini- 
mizing the cost function, and finally (v) using these new 
parameters to repeat the procedure from (i) until the pa- 
rameter values converge. Of course there is an intrinsic 
inaccuracy which, for a fixed number of bath sites can- 
not be reduced, and is due to the error in the fit of A(w) 
by a finite number of parameters. In principle, this can 
be systematically improved by increasing the number of 
bath sites. Of course, this is limited by the exponen- 
tial increase of the Hilbert space, which, in this case, is 
even faster due to the fact that the number of degrees of 
freedom of the superfermion space is twice the one of the 
fermion space, so this makes the effort more difficult than 
in ordinary ED for the same number of bath sites. One 
should, however, observe that the number of fit parame- 
ters is larger than in the case with "simple" ED without 
Markovian bath 162 . 



We have used here an effective bath containing AT 6 = 2 
sites. Still, by taking all possible parameters into ac- 
count, i.e., allowing, for example, all matrix elements of 
the r matrices to be nonzero (within the constraints im- 
posed by symmetries), this gives 8 independent fitting 
parameters. In "simple" ED, one would have only 2 in 
the particle-hole symmetric case [63] . We take model pa- 
rameters [51] v 2 = 0.1, t = 1, and the leads are fixed to 
zero temperature. In Fig. [2] we show the result of the 
fit to the bath hybridization function for two values of U 
and bias voltage $ = 2. As one can see, already with this 
small number of bath sites, the fit is quite good. More- 
over, the structure of A off is clearly non-Markovian, as ex- 
pected. For comparison, we also plot the result of the fit 
to the bath hybridization function obtained with 7V{, = 4 
bath sites. This shows a considerable improvement. The 
quality of the fit can be also inferred by directly plotting 
the local and the impurity Green's function in Fig. [3] for 
two different values of U and $ = 2. 



The bias voltage induces a steady-state current. The 
expression for the current density j (current per square 
plaquette) is obtained straightforwardly in terms of the 
layer Green's function and the v = lead Green's func- 



(a) 
0.6 



0.2 



(b) 



o U=0 
U=1 

U=2 8 
» U=3 8 
o U=4 8 
n 



8 § § § g „ 

o 

o o o o „ 



o O O o 



a o 



o o 



o O O O O o 



8S8 



0.6 ■ o o o 



0.2 



o o 







o ° 

° o • 
oo 

^ 0.4 f ° ° o o o o r 



o V =0.5 

o v 2 =0.2 
v 2 =0.1 
v 2 =0.05 



o 8 8o , 



j_£ i 







u 



Figure 4: (Color online) Scaled current density j/v 2 as a func- 
tion of (a) the bias voltage $ for different values of U and 
v 2 — 0.1, and (b) the interaction U, for different values of the 
hybridization v, and $ = 2. 



tions (see, e.g. [52]) as: 
dkii 



cIlo 



Re [G fl (k||, w ). 9i A '(k||,c) 



BZ (2tt)2 J 2tt 

+G K (k h u)gf(k h u J )] . 

This is plotted in Fig. [4] as a function of bias voltage. 
The current, as expected, decreases with increasing U for 
smaller biases. At larger $ the behavior is opposite, since 
j extends over a range of voltages which increases with 
increasing U. While at U = 0, a particle going through 
the interface conserves kii , and thus the current goes to 
zero at a bias voltage equal to the one-dimensional (z- 
direction) bandwidth, the scattering at nonzero U mixes 
k|| and thus broadens the bandwidth of possible final 
states. In Fig. [4] we plot the scaled current j/v 2 as a 
function of bias voltage for different values of U, as the 
conductance is expected to behave as v 2 in a conductor, 
while it is suppressed (cx w 4 ) in a gapped system. The 
crossing of the curves around U ~ 4 is a signal of the 
nearby equilibrium metal-insulator transition. However, 
this should be seen only as indicative, as the curves are 
taken at a relatively high bias. 

In conclusion, we have introduced a versatile method 
to deal with strongly correlated systems out of equilib- 
rium within dynamical mean-field theory. The DMFT 
self-consistent bath is approximated by an effective one 
consisting of a small number of sites coupled to a Marko- 
vian bath environment. The steady state and Green's 



5 



function of the effective system is solved by ED of the cor- 
responding Lindblad equation. The approach is particu- 
larly appropriate to deal directly with the steady state, 
without the need to consider full time evolution. Nev- 
ertheless, it should be straightforward, although compu- 
tationally more demanding, to deal with time-dependent 
problems, e.g., to describe pump-probe processes. 

The accuracy of the effective bath to reproduce the 
DMFT one obviously depends on the number of bath 
sites, which is limited by the exponential increase of the 
"super" -Hilbert space. Improvements can possibly go 
along solving the Lindblad problem in the "smaller" or- 
dinary fermion space in combination with quantum tra- 
jectory methods [65H67] . and/or by density matrix renor- 
malization group approaches [55HTD] . 

The approach illustrated here for a simple but exper- 
imentally relevant |51j model can be extended straight- 
forwardly to a number of other physically relevant sys- 
tems, including multi layer semiconducting heterostruc- 
tures, ultracold atoms and correlated coupled-cavity ar- 
rays featuring driving and dissipation, molecular con- 
tacts, and can be used to study nonequilibrium quantum 
phase transitions in these systems. Extension to a non- 
local self energy, as in cluster DMFT, or in nonequilib- 
rium variational/perturbative cluster approaches |71H73j 
is also an interesting development. 

We acknowledge illuminating discussions with S. Diehl. 
This work is partly supported by the Austrian Science 
Fund (FWF) F4103-N13, J3361-N20, and P24081-N16. 
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I. SUPPLEMENTAL MATERIAL 

In this supplementary material we present details for 
the solution of the effective bath problem described by 
the Lindblad equation 



and the evaluation of the Green's functions. 

A. Mapping onto a "super" Fock space 

Following, e.g. [56], we express the Lindblad operator 
C = Co+Ci within a "super" Hilbert space H su p = 'H®'H 
given by the tensor product of the normal % and a "tilde" 
counterpart H. Accordingly, one introduces fermionic 
operators c„, as well as their "tilde" counterparts c n . 
Within a matrix-vector notation whereby 




and similarly for c (for simplicity of notation, the spin 
index is omitted). The = (c)t are row vectors, and 
r]/2,E are the dissipation and "hopping" matrices in- 
troduced in the text. The noninteracting part can be 
expressed •">(! as 

ic = [ ct ( E - *( r i - r a )) c - (e + »(Ti - r 2 )) c 

- 2 (c T riC + c t r 2 c tT ) - 2 i trr 2 j , 

where T denotes transpose. The contribution from the 
interaction is expressed as 

iCi = Uc ^c ^c ^c Q ± — Uc ^c ^c ^c Q ±, 

B. Steady state and expectation values 

C, which plays the role of a non-hcrmitian "super" 
Hamiltonian, can be diagonalized within the many-body 
Fock-space Hsup with conventional methods. The time 
dependence of the density matrix | p ) , which is repre- 
sented as a Fock state in Hsup, is then expressed in terms 
of the (left) eigenstates ( aL | of £ with eigenvalue C a as 

(aL\p(t)) =e tc * (aL\p(0)} (5) 

Since iC is non hcrmitian, C a have real parts, which 
can be shown to be all negative (or zero) , so they produce 
decaying behavior in[5j (At least) one eigenvalue C aQ — 
0. The corresponding eigenvector | olqR ) is the steady 
state density matrix. 

Following [56], expectation values and traces are ex- 
pressed as 

{6)=tr6p = {i\6\p) = {i\dp\i) 

where 1 1 ) is the so called "left-vacuum state" [55] . No- 
tice that since = trp, one has (I \C \ p) — for any 
| p ). Thus, if the steady state is unique, one should have 



p = Cp , 



(4) 



(I\ = (a L\ 
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Conservation laws, as usual, reduce the size of the 
many-body Hilbert space one has to diagonalize. For 
additive conserved quantities, for example n a (a =t,l), 
the density matrix "state" is characterized by the con- 
straint n a — n a = 0. In general, one can define "sectors" 
with a given n„ — n^, which can be separately diagonal- 
ized. For the Green's functions, these are characterized 
by n a - n a = ±1 

C. Green's function 

To evaluate the Green's function, one needs two time 
evolutions, i.e. expectation values of the form [55] 

G(B(t),A(t')) = truB(t)A(t')pu = trB A t ,{t - t') (6) 

for t > t', where pu is the density matrix of the universe 
(in contrast to p, which is the reduced density matrix of 
the "system"), and try = tr ® trg is the trace over the 
universe degrees of freedom, whereby tr is the one over 
the system, and trg the one over the environment. Here, 

A t ,(s)=tr E e- isH "Ap u (t')e isH " 

and Hjj is the total Hamiltonian of the universe. 

The Quantum regression theorem |58j states that under 
the same assumptions for which Q holds, one has (for 
s > 0) 

±M*)=£M*)- (7) 

as 

In terms of eigenvalues and eigenvectors, this gives 



For example, the "greater" Green's function 
iG >ntTn {t,t') = tvuc^cKt^pu =(I\ c4M) I 1 ) 

is obtained from ^ with B — c n , and A — c\ n . Inserting 
the identity ^ Q | aR ) ( aL | and taking the steady state 
| pit! — +oo) ) = | aoR ) this gives 

iG> n , m (t' + s, t) = 

= e (*-*')^« ( ao L | c„ | aR ) ( aL | c} n \ a R ) 

a 

for t > t' . In a similar way, one obtains the result for 
t < t' and for the "lesser" Green's function. These are 
combined (and Fourier transformed), to obtain the re- 
tarded (G r (oj)), advanced G a (uj) 1 and Keldysh G K [uj) 
components of the Green's function of the effective bath 
problem evaluated at the impurity site {n — m — 0). 

The self-energy is extracted by using Dyson's equation 
in the 2x2 Keldysh space: 

S(w) = G^lu)- 1 -Giuj)- 1 . 
Using a modified version of the "Langreth rules" [52] , one 
can identify the separate components 

= 1/G » - 1/G» , 

-£ a (lo) = Z r (uj)*, and 

£ A » = -G«(u)/\G R (u)\ 2 + G K (u>)/\G R (u)\ 2 . 



(aL\A t ,{s)) = e sC ~ (aL\A t> {0)) 



